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ABSTRACT 

We use the results of large-scale simulations of reionization to explore methods for charac- 
terizing the topology and sizes of HII regions during reionization. We use four independent 
methods for characterizing the sizes of ionized regions. Three of them give us a full size dis- 
tribution: the friends-of-friends (FOF) method, the spherical average method (SPA) and the 
power spectrum (PS) of the ionized fraction. These latter three methods are complementary: 
While the FOF method captures the size distribution of the small scale H II regions, which 
contribute only a small amount to the total ionization fraction, the spherical average method 
provides a smoothed measure for the average size of the H II regions constituting the main 
contribution to the ionized fraction, and the power spectrum does the same while retaining 
more details on the size distribution. Our fourth method for characterizing the sizes of the H 
II regions is the average size which results if we divide the total volume of the H II regions 
by their total surface area, (i.e. 3V/A), computed in terms of the ratio of the corresponding 
Minkowski functionals of the ionized fraction field. To characterize the topology of the ion- 
ized regions, we calculate the evolution of the Euler Characteristic. We find that the evolution 
of the topology during the first half of reionization is consistent with inside-out reionization 
of a Gaussian density field. We use these techniques to investigate the dependence of size and 
topology on some basic source properties, such as the halo mass-to-light ratio, susceptibility 
of haloes to negative feedback from reionization, and the minimum halo mass for sources to 
form. We find that suppression of ionizing sources within ionized regions slows the growth 
of H II regions, and also changes their size distribution. Additionally, the topology of sim- 
ulations including suppression is more complex, as indicated by the evolution of the Euler 
characteristic of the ionized regions. We find density and ionized fraction to be correlated on 
large scales, in agreement with the inside-out picture of reionization. 

Key words: H II regions-ISM: bubbles-ISM: galaxies: high-redshift-galaxies:formation- 
intergalactic medium-cosmology:theory 



1 INTRODUCTION 

The Cosmic Microwave Background (CMB) discovered in 1965 
vyas evidence that the hot big bang universe cooled and recombined 
jPenzias & Wilsonll965l) . That same year, however, the intergalac- 
tic medium (IGM) at z = 2 was found to be largely devoid of neutral 
hydrogen atoms, when astronomers failed to detect their Lyman al- 
pha resonant scattering in the spectra of the first quasars discovered 
with hi gh enough redshift to m a ke the tran sition visible from the 
ground dOunn & Petersonlll965l : lok3ll966h . This was soon inter- 
preted to mean that, unless the IGM were many orders of magni- 
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tude less dense than the average density of a critical universe, most 
of the hydrogen atoms there must have been reionized sometime 
between z = 1000 and 2 = 2. 

Although we have since then learned much more about both 
the CMB and the HI absorption towards high redshift QSOs, cur- 
rently it is still those two observables which constrain the epoch of 
reionization (EoR). The results from the WMAP measurements of 
the CMB have constrained the optical depth due to electron scatter- 
ing, Tcs, to 0.088 ±0.015, implying that an i nstantaneous reioniza - 
tion would have happened at 2: — 10.4±1.2 jKomatsu et alj201oh . 
The QSO spectra obtained within the Sloan Digital Sky Survey 
(SDSS) indicate a low, but rapidly rising neut ral fraction around 
redshift 6 tFan et al.i,2006 : IWillott et alj|2007h . The combination 
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of these two measurements suggests that the epoch of reionization 
extended over several redshift units. 

However, a series of measurements are being performed or 
prepared that are expected to give completely new constraints on 
this early epoch of galaxy formation. Radio telescopes capable 
of measuring at low frequencies (GMRtB 21CMA0, LOFAF0, 
MW4E PAPEfQ) should be able to detect the signature of red- 
shifted 21cm radiation from neutral hydrogen during the EoR. 
These measurements are challenging due to the presence of the 
strong foreground emission of, mostly, our own Milky Way as well 
as ionospheric distortions. If successful, these experiments should 
produce a detection before the 50th anniversary of the discovery of 
the CMB and the Gunn-Peterson effect. 

In preparation for these 21cm observations, many groups 
have been nu merically s imulating the reionization process on 
large scales (e.g. llliev et a l. 2006a; McQuinn et al. 2007; Shin et al.l 
I2OO8L to name a small selection); see iTrac & GnedinI (|2009j) for 
a review on simulations of reionization. Semi-numerical models 
have also been developed, first 
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by IZahn et alJ (l2007t). later im- 
rlanetto"2007';'santos et al."2008'; 



proved by others ( Mesinger & Furlanetto 
Alvarez et al. 2009; Choudhurv et al. 2009; Mesinger et al. 2010; 
Alvarez & Abelll2010h . What both of these types of calculations 



give us is the evolution of the ionized fraction of intergalactic hy- 
drogen in the Universe, a;(r, t). 

The simulation results show a great amount of complexity in 
x{r, t). As the sources of reionization are likely to be clustered in 
space, individual H II regions typically contain many sources (e.g., 
iFurlanetto et al.l2004l ; riliev et al.l2005l.l2006bl) and obtain complex 
shapes in 3D space. Rare sources are more biased than more abun- 
dant ones, and it is expected that the level of bias will larg ely deter- 
mine the characteristic scale o f the reionization process ( Iliev et al.l 
l2006bl ; IFurlanetto et alj|2006h . Accurate theoretical predictions for 
the morphology and size of H II regions depend upon an under- 
standing of the the abundance and clustering of the ionizing sources 
themselves, in addition to the underlying inhomogeneous density 
field. Quantitative analysis of the distribution of ionized material 
during the EoR is thus not a trivial matter and likely several differ- 
ent approaches have to be combined. The main aim of this paper is 
to describe and evaluate different methods for analyzing the prop- 
erties of the ionization fraction field x{r, t) and its evolution. 

The planned observations of the 21cm line of neutral hydrogen 
are expected to constrain the ionization fraction field x{y, t) statis- 
tically as the power spectrum of neutral hydrogen fluctuations is the 
most directly observed quantity via 21-cm ra di o observ atio ns (e.g . 
IZaldarriaga et alj|2004 fMellema et alj|2006l ; fHarker et al]|201ol) . 
Future observations, for example with the SKA (Square Kilometer 
Array) may have enough sensitivity to actually image the redshifted 
21cm signal as a function of frequency and thus reveal the spatial 
structure of the ionization fraction field x{r, t). 

Ultimately, we are less interested in the function x{r, t) itself 
but more in "why it is like it is", i.e. in the properties of the sources 
and sinks of reionization. For this, one has to find out how the sta- 
tistical properties of ionization fraction field depend on different 
source and sink properties. Simulations of reionization can thus be 
said not to aim at reproducing the actual x(r, t), but rather at show- 
ing the same statistical behavior as the real epoch of reionization: 
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Figure 1. Evolution of the effective efficiency factor g-y for the simulation 
with imposed photon production history (53Mpc_uvS_le9) as a function of 
time and global ionization fraction. 



on the one hand, because only statistical quantities can be mea- 
sured (as mentioned above), on the other hand because the input of 
simulations is only statistically comparable to the real conditions. 

This paper has two parts. In the first part we investigate the 
usefulness, in terms of characterizing size distributions of ionized 
regions, of different kinds of statistics (for example the power spec- 
trum of x(r, t)) of a simulated ionization fraction field. In the sec- 
ond part we employ these statistics to investigate the effect of dif- 
ferent source properties. This is useful to draw conclusions on the 
sources, once statistical properties of the real x{y, t) of reionization 
can be measured. 

We focus on the early and intermediate stages of reionization, 
when the morphology of H II regions is most well defined and the 
photon mean-free-path is determined by the patchiness of the reion- 
ization process itself. At the latest stages of the reionization pro- 
cess, after overlap, fluctuations in the UV background are expected 
to be sensitive to the small fraction of gas whi ch is left neutral in the 



20001; iGnedin & Faij|2006l ; I Alvarez & Ab"eill2010l ; iProchaska et al 
120091) . We limit ourselves to analyzing the ionization fraction fields 
x{r, t) from the simulations, not on producing the observable quan- 
tities. This is the necessary first step before proceeding to evaluate 
whether different scenarios can be observationally distinguished. 
The observables will be discussed in a follow-up paper (Iliev et al., 
in prep.). 

In terms of sections, the paper is organized as follows: §2 in- 
troduces the simulations included in this study. §3 introduces the 
analysis methods used to investigate these simulations. In §4 we 
test the effect of numerical parameters on the statistics of x{r, t). 
In §5, we test the effect of source properties on the statistics of 
x{v, t). We end with our conclusions in §6. 



2 SIMULATIONS 

Our simulation m ethodolo gy has been previously described in de- 
tail ( Illievetal.l20 06b; Mel lema et alj2006l : llliev et alj2007l) . Here, 
we will briefly summarize the underlying N-body simulations that 
were performed and the set of radiative transfer simulations that 
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Table 1. Simulation parameters, volumes derived from this and global (mass averaged) reionization history results for simulations with WMAP5 cosmology 
parameters. The box sizes can be directly inferred from the simulation names. For all simulations, the mesh consists of 256^ cells. The ionization time step 
for all simulations is At^ = 5.75 X 10^ yr. g-^ is the efficiency parameter as explained in the text; the old efficiency paramter is given in brackets; Vmin 
is the comoving volume of the minimum size H 11 region, ionized by the least efficient source during one time step, assuming the density to be the average 
density of the universe; Tcs is the electron scattering optical depth calculated for each simulation. 

53Mpc.g8.7.130S 163Mpc.g8.7.130S 53Mpc.g8.7.130 53Mpc.gl.7.8.7S 53Mpc.g0.4.5.3 53Mpc.uvS.le9 53Mpc.gl0.4.0 



miss g-r (h) 8.7 (10) 

mrssff7 (/7) 130(150) 

suppression yes 



8.7 (10) 

130(150) 

yes 



8.7(10) 
130(150) 



1.7(2) 

8.7(10) 

yes 



0.4 (0.4) variable -^Fig[T] 10.4 (12) 

5.3 (6) (0) (0) 

no n/a n/a 



Vccii/Mpc-=' 

Vmin/Mpc3 



0.0088 
0.1361 



0.2575 
0.1361 



0.0088 
0.1361 



0.0088 
0.0136 



0.0088 
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Figure 2. Ionization maps (ionization fraction according to colour bar) of all included simulations at 30% global ionized fraction: from top left to bottom light: 
53Mpc.gl0.4.0, 53Mpc.uvS.le9, 53Mpc.g8.7.130. 163Mpc.g8.7.130S, 53Mpc.g0.4.5.3. 53Mpc.g8.7.130S and 53Mpc.gl.7.8.7S. Each panel is for a slice 
which is one cell thick (~ 0.64 Mpc and ~ 0.21 Mpc, respectively for the 163 Mpc and 53 Mpc simulations). 
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we analyze. On the scales of interest to us here, the intergalactic 
medium and dark matter followed each other as cosmic structure 
arose in the ACDM universe during this epoch. Even during reion- 
ization, since ionization fronts which reionized the IGM moved 
supersonically, the back reaction of the gas due to mass motions 
related to pressure force s can be neglected to first approximation 
dShapiro & Gir oux 1987), and, hence, the radiative transfer can be 
done as a post-processing of the N-body density field. 



2.1 N-body simulations 

As a basis for our radiative transfer calculations, we begin with 
the time-dependent density field extracted from N-body simula- 
tions of structure formation. We use the CubeP'^M code w hich 
was developed f rom the PMFAST code ilVIerz et all |2005|) . see 
llliev et al.l ( |200^ for a short description of the CubeP^M code. 
It uses particle-particle interactions at sub grid distances and a 
particle-mesh method for larger distances. Here, we use the re- 
sults of two simulations, performed with CubeP^M, one for a vol- 
ume of 163 Mpc on a side, another with 53 Mpc. The former has 
3072'^ particles and mesh size of 6144^ cells, while the latter has 
1024"^ particles and 2048^ cells, which imply particle masses of 
5.5 X 10'' M0 and 5.1 x lO'' Mq, respectively. These parame- 
ters guarantee a minimum resolved halo mass of 10* M0 which IS 
approximately the minimum mass of halos able to cool by atomic 
hydrogen cooling. The cosmological parameters used were for a 
flat ACDM universe with fi™ = 0.27, Sit = 0.044, h = 0.7, 
n = 0.96, and era = 0.8, based on the five year WMAP results 
jKomatsu et"ai]|2009l) . 



2.2 Radiative transfer runs 

Table [T] gives an overview of the seven different radiative transfer 
runs that we analyze in this paper. These are a sub-set of a larger 
suite of simulations, to be presented in a follow-up paper (Iliev et 
al., in prep.). This sub-set was chosen as the minimum one needed 
to illustrate the points we want to make in this work. All the radia- 
tive transfer simulatio ns were performed using the C^-Ray method 
jMellema et al.l2006l) on a uniform rectilinear grid containing 256 
grid cells. The density is assigned to the mesh by smoothing the 
dark matter particle distribution from the underlying N-body sim- 
ulation using an SPH kernel function: each N-body particle is as- 
signed a compact, spherical smoothing kernel whose width is ad- 
justed so as to encompass its 32 nearest-neighbors. Particle mass is 
then assigned to the cells of our radiative transfer grid by integrat- 
ing each kernel function over the volume of each cell it overlaps. 
To convert what is a the IGM dark matter density into a baryon 
density, we assume that in the IGM the gas distribution follows the 
dark matter. This is valid on the scales of the radiative transfer cells 
(0.2 or 0.6 comoving Mpc) as at the mean density of the IGM they 
are much larger than the local Jeans length. 

There are physical effects below our resolution limit which 
influence reionization. These are small scale density variations 
('clumping') and the presence of unresolved absorbers (such as 
mini-halos and the structures of unknown origin which are ob- 
served as Lyman Limit Systems at lower redshifts). All of these will 
slow down the reionization process as they increase the number of 
photons absorbed. In the simul ations prese nted here we d o not con- 
sider t h ese effects, but see e.g.lCiardi et al. (2006): McO uinn et al.l 
( l2007h : lAlvarez & Abeil j2010t) : ICrociani et ali 12010) for studies 
about the effect of different types of unresolved absorbers and e.g. 



llliev et al.l ( l2006bh ; iMcOuinn et al] feOOTh : llUev et alj ( l2008h for 
compartative studies of the effect of clumping. 

Simulations are labeled with the parameter g-y, which is an 
efficiency factor for the ionizing photon production of halos per 
source halo baryon per unit time. Each halo of mass M is assigned 
a luminosity 



N, 



(1) 



where N-y is the number of ionizing photons emitted per Myr, M is 
the halo mass, and rrip is the proton mass. Halos are assigned dif- 
ferent luminosities according to whether their mass is above ("large 
sources") or below ("small sources' ') 10*'Mq (but above l(f Mq). 
For example, 53Mpc_g8.7_130S indicates that large sources have 
an efficiency = 8.7, while small sources have an efficiency 
g~f — 130, and the symbol "S" means that the small sources are 
suppressed in regions where the ionization fraction is higher than 
10%. 

In previous simulations performed with C'^-Ray, the source 
efficiencies were characterized by f-, , the number of ionizing pho- 
tons released per source halo baryon per star-forming episode (i.e. 
per simulation time-step for updating the source halo catalogue 
from the N-body results). The relation between f-y and g.y is given 
by 



(2) 



where At is the time between two snap shots from the N-body 
simulation. For example 'liievetai] IIOO^ considered a simulation 
called f250. In the new naming scheme this would be called gl 12. 
The reason for switching to a new naming scheme is that the pre- 
vious scheme hid the dependence on the size of the time step At, 
since ionizing photons were released over a time At per baryon 
for all the baryons inside source halos when that step began. This 
made it more difficult to compare simulations involving different 
time-steps, since the results depend on BOTH f-, AND At, while 
the instantaneous luminosities of source halos depend only upon 
their ratio f^/At, not alone. 

The suite of simulations presented in Table 1 allows us 
to see how the morphology and characteristic scales of reion- 
ization depend upon various important numerical and physical 
parameters which are not yet well understood. The simulation 
53Mpc_g8.7_130S is our standard case for this paper. We refer to 
this as the fiducial simulation. It produces an electron scattering op- 
tical depth consistent with the l-cr ra nge allowed by the se ven year 
WMAP results, Tes = 0.088 ±0.015 jKomatsu et alj20IOh . To test 
the effect of weaker sources, and thus more extended reionization 
we also present 53Mpc_gl.7_8.7, which ends considerably later and 
has an optical depth consistent with the seven year WMAP results 
when considering the 2-a range (and assuming a gaussian error 
distribution) for Tcs- These two simulations are used to introduce 
the different analysis methods in section[3] 

One of the physical effects which may be present during 
reionization and which we study in this paper, is source suppres- 
sion due to Jeans mass filtering, in which ionizing radiation from 
sources hosted by halos with a mass below some threshold is sup- 
pressed when the h alos are located within ionized regions (e.g. 
IShapiroetai]|l994l) . This concept was introduced in our simu- 



Iliev et al. I ( l2007h . By comparing, for example 



lation models in 
53Mpc.g8.7.130S to 53Mpc_g8.7_130 (the latter without source 
suppression), it is possible to isolate the effects due solely to source 
suppression. However, the simulation with no suppression will end 
at a much higher redshift and therefore the halo populations are not 
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comparable at corresponding stages of reionization (e.g., at 2:50%). 
Hence, we also include the simulation 53Mpc_g0.4_5.3 which does 
not have suppression, but which due to the weaker source luminosi- 
ties, ends approximately at the same time as our fiducial simulation 
53Mpc_g8.7_130S. This way the different reionization stages (ex- 
cept the earliest ones) occur at similar times, and thus these two 
simulations have similar halo populations at the different stages of 
reionization. 

For the 54 Mpc simulation volume at z ~ 13.6, there 
are roughly 330 cells containing source halos more massive than 
IO^Mq (this corresponds to roughly 2 x 10"^ % of all cells). 
Additionally, there are about 38 000 cells (2.5 x 10"^%) contain- 
ing low mass source halos between lO^M© and IO^Mq. However, 
for the fiducial simulation for example, roughly 88% of these low 
mass source halos are suppressed. For the large simulation volume, 
these numbers (also at roughly 10% global ionization fraction, i.e. 
z ~ 13.2) are: 12 000 cells containing massive halos (7.8 x 10~^% 
of all cells) and 750 000 cells (4.8% of all cells) containing low 
mass halos, of which roughly 86% are suppressed. At overlap, the 
small (large) simulation volume has about 17000 (440 000) cells 
containing massive source halos, which corresponds to 0.1 % and 
2.8 % of all cells, respecively. The number for low mass halos are 
280 000 (3 600 000) cells, or 1.8 % and 13 % of all cells, almost all 
are completely suppressed. 

Throughout this study, we will make comparisons like indi- 
cated above for the case of source suppression, in order to see how 
physical effects manifest themselves and to find which quantita- 
tive measurements best discriminate among different reionization 
scenarios. Instead of comparing the simulations at equal redshifts, 
we do the comparison at equal mass averaged (global) ionization 
fraction, (a;). Table 1 lists the redshifts at which the global ioniza- 
tion fraction (i.e. mass-weighted average, unless otherwise stated) 
for each simulation are (a;) ~ 0.1, 0.3, 0.5, 0.7, 0.99. The epoch 
at which the H II regions globally "overlap", Zov, will, by conven- 
tion, be taken here to be the redshift at which (2;) = 0.99, although 
the value of 2:0V which results is not very sensitive to this particular 
choice as long as (a;) is close to unity. 

Beside the above mentioned simulations, there are three 
more simulations included in this study: the simulation labeled 
53Mpc_uvS_le9 has the same (imposed) global photon production 
history as our fiducial simulation, but only halos more massive than 
are allowed to host luminous sources. This results in a g-, 
that is time dependent. Its evolution is plotted in Fig.[T]as a function 
of time and mass averaged ionization fraction (a;) . 

As a second simulation with only high mass sources, we in- 
clude simulation 53Mpc_gI0.4_0. Unlike 53Mpc_uvS_Ie9 it has a 
constant mass to light ratio which is chosen so that reionization 
ends roughly at the same time as in our fiducial simulation. This 
yields later 210% — 270% and thus a lower value of Tcs. Note that 
the efficiency of sources in high mass halos had to be boosted only 
by a factor 1.2. This means that in our fiducial simulation, sources 
in low mass halos only contribute about 17% to the ionizing photon 
budget. 

Simulation I63Mpc_g8.7_I30S has the same physical param- 
eters and the same mass resolution and, hence, halo mass range as 
our fiducial simulation, but the simulation box volume is about 30 
times bigger. Therefore, it is capable of catching structure on larger 
scales. On the other hand, the resolution in the radiative transfer 
simulation is worse than in the small box simulation since the num- 
ber of cells in both simulations is 256 per side. This simulation is 
included to check for cosmic-variance effects and to test the effect 
of resolution on our investigation methods. 



Snapshots of the simulations at the 30% global (mass aver- 
aged) ionized fraction are shown in Figure (2] The slices are to the 
same comoving physical scale to make it more easy to see the mor- 
phological and topological differences between the models which 
we will discuss in detail below. 

Table [T] also lists values for the smallest possible H II region 
Vmin which could be formed during a single radiative transfer time 
step Ati if the surrounding IGM has the average density of the Uni- 
verse and recombinations can be neglected. This number is likely 
to be an overestimate as recombinations and density peaks will re- 
duce it. However, it is a useful number to compare the resolution of 
various radiative transfer simulations with. The number of emitted 
ionizing photons from the smallest halos of mass A/min is given by 

N ■ - N At - gT Atj 
10 fiiTipilo Myr 

which then gives a minimum volume of 

riH fiopcrit.o lOMyr 
- 0.0016Mpc3(^),.. (4) 

where we have used the radiative transfer time step Ati — 5.75 x 
10*' yr and nn = 0,b is the hydrogen number density. Vmin 

can be compared to the cell-sizes of the simulations which are also 
listed in Table[T] 



3 INTRODUCING THE ANALYSIS METHODS 

In this section we introduce our analysis methods by means of two 
simulations that differ only in the mass-to-light ratio of the ha- 
los, 53Mpc.g8.7_I30S and 53Mpc_gI.7_8.7S. Whenever we refer 
to "the simulations" in this section, we mean these two. The focus 
in this section is on the ability of our analysis methods to discrim- 
inate between the two simulations. The results of all methods to- 
gether can be seen as a characterization of the morphology of the 
ionization fraction. 

3.1 Size distribution 

One of the most basic measures of reionization is the size distri- 
bution of H II regions. However, as will become clear below, "size 
of an H II region" is a quantity which can be defined in differ- 
ent ways. Under the assumption that most of the volume is either 
highly-ionized or highly-neutral, H II regions can be considered to 
be topologically connected volumes of spac e. We previously used 
a friends-of-friends (FOF) method jlliev et al. 2006b) to identify 
such regions, using the condition x > 0.5 for a cell to be con- 
sidered ionized. In co ntrast to thi s meas ure of the volume of con- 
nected ionized space, IZahn et al. used a different method, 
introduced as "the bubble probability distribution". For reasons we 
explain below, we refer to this method as the "spherical average" 
method. We now describe these two methods in more detail. 

3.1.1 Friends-of-friends method 

Our first method for identifying the size distribution of H II regions 
relies on a literal definition of "H II region": a connected region in 
which hydrogen is mostly ionized. For grid data, the obvious way 
to identify such a connected region is to use a "friends-of-friends" 
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(FOF) approach, in which two neighbouring cells are considered 
friends if they both fulfill the same condition. Cells are grouped 
into distinct regions according to whether they are linked together 
in an extended network of mutual friends. The algorithm we use to 
group cells togeth er is the equivalence class method, described in 
IPress et aLl l ll992h . Unless otherwise specified, we use x > 0.5 for 
a cell to be considered ionized, and x < 0.5 for a cell to be consid- 
ered neutral, so that every point in the simulation box is either in an 
H I or a n H II region. Our method was first described in llliev et al.l 
(l2006tj). In contrast to all other size measures presented below, the 
FOF does not care about how contorted an H II region is. Therefore, 
the sizes of H II regions found by the FOF method and the sizes of 
H II regions found by other methods which will be introduces be- 
low, give complementary information about the morphology of the 
ionization fraction field. 

The FOF method has been use d extensively for h alo finding in 
cosmological N-body simulations ( Davis et al.l l 198? ). Our imple- 
mentation is more straightforward, since each cell always has only 
6 direct neighbours, the identities of which are known in advance, 
as opposed to particle data, in which it is necessary to perform 
costly searches to identify the groups. Another significant differ- 
ence between the two methods is the role played by free param- 
eters. In the halo finding FOF method, the free parameter is the 
linking length, which is the distance within which two particles are 
considered to be friends. In the region finding method, the free pa- 
rameter is the threshold, a;th, for a cell to be considered ionized or 
neutral. 

As seen in Fig. [3] we test the effect of varying xth (using 
our fiducial simulation 53Mpc_g8.7_I30S at 50% global ionization 
fraction), where we used three different values for a::th ~ 0.1, 0.5 
and 0.9. The thin dot-dashed black line shows at every volume V, 
how much a single region with this volume would contribute. This 
means that each bin cannot be filled less (i.e. or the bin must be 
empty) than to the point where this line crosses the lower limit of 
each volume bin. 

As can be inferred from Table [T] the minimum volume ion- 
ized by a source in a single radiative transfer time step for our 
fiducial simulation is about an order of magnitude bigger than the 
cell size. This means that at the end of the radiative transfer time 
step at which the source turned on, the cell hosting the halo will 
be completely ionized and surrounded by partly ionized cells. If 
regions are connected through such partly-ionized border cells, it 
strongly depends on xth if the regions count as two disconnected 
or as one connected region. In every analysis method that depends 
on a threshold value this effect is bigger if the partly ionized bor- 
ders of ionization regions are comparable to the regions' size. Nev- 
ertheless, as can be seen in Fig. |3] the qualitative picture remains 
unchanged, with a few few-cell regions, a substantial contribution 
of regions of intermediate size, and the main contribution coming 
from a single large region comparable in size to the simulation box. 
The absence of single-cell regions for a;th = 0.1 can be explained 
by looking at the minimum number of photons from a single source 
released during one time step. Sufficient photons are produced to 
ionize the cells surrounding the source cell more than 10%, even 
if the density of the cells is nine times the average density of the 
universe, see Eq.|4] 

Below, we convert the volume bins into equivalent radius bins, 
-Rcqui ,FOF = [3/{4:7t)V]''^^'^^ which corresponds to the radius 
^equi ,FOF of a spherical region with the same volume. We want 
to stress that this does not mean that the H II regions are spheri- 
cal. We convert to an equivalent radius only to allow a more direct 
comparison to the spherical average method and the power spec- 
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Figure 3. Effect of varying the thre.shold Xth for the FOF method for the 
53Mpc_g8.7_130 simulation at z = 9.7, when the global ionized fraction 
(x) is about 50 percent. The results for three different threshold limits are 
shown, as indicated in the figure. The thin dot-dashed black line shows at 
every volume V, how much a single region with this volume would con- 
tribute. Indicated by " on the abscissa ai'e the cell- and box- volume. 

trum which are described in the following sections. We normalize 
to the total volume and not to ionized volume, to allow for a more 
direct comparison to the power spectrum which is normalized in 
the same way. 

In order to show the time-evolution of the FOF size distribu- 
tion of ionized regions for a simulation in a single plot, we choose 
a fixed threshold value (xth ~ 0.5) and color code the contribution 
Vdp/dV. The color coding makes it possible to show a histogram 
(like Fig. [3]( in a single line or column. Each individual column 
of Fig. |4] (left panel) is a histogram as Fig. [3] at a different global 
ionization fraction. Fig.|4]thus shows the evolution of the size dis- 
tribution with global ionization fraction. The box- and the cell size 
are marked with > on the ordinate. By construction the FOF method 
will not result in sizes larger than the former and smaller than the 
latter. As we will see later, H II regions start to merge very early 
on in the course of reionization which results in shapes far from 
spherical and a complex topology already at global ionization frac- 
tions (x) ~ 0.2. The concept of distinct H II regions and their sizes 
quickly becomes meaningless. Therefore, all size distribution esti- 
mates are only shown up to a global ionization fraction (x) ~ 0.6. 

Three things catch the eye when analyzing the size evolution 
in FigH 

(1) Already at (x) ~ 0.15, the distribution for both simula- 
tions is not continuous, but shows a gap. Most of the ionized vol- 
ume is contained within one region of a size falling into a size-bin 
which is separated from the rest. This is an inherent property of the 
FOF method, where regions are grouped together as soon as they 
touch and local H II region percolation occurs quite early in the evo- 
lution. If the H II regions reach a certain size, which depends on the 
clustering of the sources and on their efficiency, they will percolate 
and form bigger H II regions. As those smaller H II regions grow 
and merge into the larger one, both their numbers and the fraction 
of the ionized volume that they occupy decrease. A doubling of 
the volume (merging of two bubbles with the same size) thus cor- 
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Figure 4. Size distributions using the friends-of-friends method for simu- 
lation 53Mpc.g8.7.130S (53Mpc.gl.7.8.7S) in the left (right) panel as a 
function of global ionization fraction {x). Vdp/dV is colour coded ac- 
cording to the colour bar. The equivalent radii corresponding to the cell and 
box- volume of the 53 Mpc box, iJoqui,FOF = 0.13 Mpc and 32.8 Mpc 
are marked by >, respectively. Additionally, the cell size for the 163 Mpc 
simulations, Roqui.FOF = 0.4 Mpc see section|4] is indicated by -. 

responds to a jump over one effective radius bin. The largest H II 
region grows through mergers with smaller H II regions, as well as 
due to sources within the region, and approaches the box-size by 
(a;) ~ 0.6. 

(2) The distribution of the smaller scales (i.e. everything 
except the one big region) is much flatter for the simulation 
53Mpc_gl.7_8.7S. This simulation has more single- or few-cell- 
sized regions than our fiducial simulation 53_Mpc_g8.7_130S. 

(3) While the size-bin which contributes most to the global 
ionization fraction is increasing with (x), the shape of the distri- 
bution of the rest does not change much. However, its total contri- 
bution to the global ionization fraction decreases with (x): At the 
same rate that small H II regions grow bigger and merge into even 
bigger H II regions, new small ones are "born" but contribute in the 
course of reionization less and less to the global ionized-fraction. 

The FOF-method applied to an ionization field in a finite sim- 
ulation box can only sample the true underlying size distribution 
function (as it would be in an infinite simulation box) up to a lower 
limit. This is indicated by the thin dot-dashed black line in Fig[3] 
This lower limit depends on the size of the simulation box. The gap 
mentioned above will be smaller the better the sampling, that is the 
bigger the simulation box. 

3.1.2 Spherical Average method 

The sp herical average method (SPA) was described bv lZahn et al.l 
( 12007 ): it can be easily used for comparisons with analytical mod- 
els. It is based on constructing spheres around every cell in the com- 
putational volume, averaging the ionization fraction inside these 
spheres and finding the largest such sphere for which the average 
ionization fraction is greater than a certain threshold xth- We chose 
Xth = 0-9. Because of this, we call it the spherical average method 
(SPA). It yields a smoother distribution of H II region sizes than the 
one obtained by the FOP method. It does not measure the size of a 
connected ionized space but instead it is a measure of the scales of 
spherical bubbles which would cover the ionized space. 

Motivated by the analysis given in Appendix lAl we multiply 
the radius found by the spherical average method by a factor of 4, 
s = 4 X i?. We call this the scale of the spherical average method. 
In Fig.[5]we plot the spherical average distribution in the same way 
as the FOP distribution in Fig. |4] RdP/dR is normalized to the 
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Figure 5. Size distributions using the spherical average method as explained 
in the text, log^Q RdP/dR colour coded according to colour bar; left 
(right) panel simulation 53Mpc.g8.7.130S (53Mpc.gl.7.8.7S) as a func- 
tion of global ionization fraction (x). The peaks of the distribution at every 
global ionization fraction are indicated by the dashed line. The black solid 
line shows the 3VA size measure, as explained in section [3. 1.31 

whole volume. We define i?max ((a;)) to be the position of the max- 
imum of RdP/dR at every (x). We see that i?max is increasing 
with (x) which means an increase in the average bubble size with 
global ionization fraction, as expected. Although initially smaller, 
it can be seen that the average scale of simulation 53Mpc_g 1 .7_8.7S 
is growing faster with respect to (a;) than the average scale of sim- 
ulation 53Mpc_g8.7_130S. Prom Fig.[5]it can be further seen that 
the distribution of the model with the smaller efficiencies shows 
initially (at (x) < 0.1) a wider distribution of bubble sizes with a 
peak at smaller scales: The smallest H II regions are smaller than in 
the fiducial model, but the contribution from small regions to total 
ionized fraction is smaller. Therefore the bigger H II regions have 
to grow to a bigger size in the 53Mpc_gl.7_8.7 model to reach an 
ionization fraction of 10%. At higher average ionization fractions, 
the peak of the distribution is slightly shifted towards bigger scales 
with respect to the 53Mpc_g8.7_130S model. 

As we have already noted, much of the noticeable difference in 
the POP curves comes from a very small fraction of the volume, and 
a log scale is required to see such differences in the POP plots. The 
spherical average distributions are much smoother, and therefore 
offer a less detailed, more global picture of the spatial structure of 
the ionized regions. However, the spherical average shows much 
more clearly the difference in size of the large scale H II regions 
between the two simulations. 

Another method which is similar to the SPA is a method used 

in|Mesinger & P urlanetto (2007). We examine this method in Ap- 
pendix |B] We find that it yields qualitatively the same results as the 
SPA and PS but its use in simulations with a continuous distribution 
of ionization fractions (i.e. not a binary field) creates complications. 

3.1.3 3V/A method 

As another estimate of the scale of bubbles, one could use the ratio 
of the total volume of the H II regions and their total surface area: 

3 X ^ Volume/ ^ Area. (5) 

H II regions H II regions 

The volume V and the surface area A were calculated from the 
zeroth and first minkowski functionals, respectively: V — Vo and 
A — 6 X Vi. For a distribution of disconnected spherical bub- 
bles, 3V/A is the surface weighed average radius. For three dimen- 
sional bodies, one could say that 3V/A is proportional to the sur- 
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Figure 6. Power spectra of ionized fraction as a function of global ioniza- 
tion fraction, A^^(fc) colour coded according to the colour bar; left (right) 
panel for the fiducial simulation (53Mpc_gl.7_8.7S). The dashed line indi- 
cates the peak of <5^^(fc) as a function of {x). Note that A^(fc) is plotted 
over a multiple of inverse k, for details see text. 

face weighted average of the smallest scale of each object. If the 
dominant structures are 2-dimensional, i.e. disc-like, then 3V/A is 
three times bigger than the surface weighted average disc-height. If 
the dominant structure is 1 -dimensional, i.e. bar-like, 3V/A is 1.5 
times the surface weighted average bar-radius. In any case, in terms 
of surface weighted averages, it is an overestimate of the minimum 
scales. As in the case of the spherical average, 3V/A does not mea- 
sure the sizes of topologically connected ionized volumes. 

It can be seen (Fig.O that initially, for both simulations, the 
3V/A scale agrees with the scale of the maxima of the spherical 
average. At about 20% global ionization fraction, the scale of the 
maxima of the spherical average for the 53Mpc_gl.7_8.7S simu- 
lation is greater than the scale estimated by 3V/A, indicating that 
while most of the volume is contained in larger scale bubbles, most 
of the surface comes from small scale bubbles: many very small 
scale structures, few very large scale structures. For the fiducial 
simulation, the offset between the scale of the maximum of the 
spherical average and the scale estimated by 3V/A is smaller, indi- 
cating that the same bubbles which contribute substantially to the 
total volume, contribute substantially to the total surface area. 

3.2 Power spectra 

The typical reionization scales can be further characterized by 
the appearance of features in the power spectra of the density 
and ionized fraction fields, Pgs, P^s, and Pxx, where (SkS^.,) = 
(27r)3j3(k - k')PM(fc), (Wk'> = 5'(k - k')(27r)^P.5(fc), 
and (5xk<5*k') = '5^(k - k')(27r)^P^:^(fe). Here, S is the over- 
density of matter, while 5x = x — Xv, where Xv is the volume- 
weighted global average ionized fraction. Note that we do not nor- 
malize 5x by x-v When plotting the actual power spectrum, we use 
the dimensionless power per logarithmic interval in wavenumber, 

A^{k) = k^P{k)/{2TT^). 

The power spectrum of the ionization fraction field is one 
componen t in the expansion of the neutral hydrogen power spec- 
trum (e.g. lPurlanetto et aill2006h . As pointed out in section 1, fu- 
ture 21-cm radio observations aim at observing this quantity via 
the 21cm line of neutral hydrogen. Therefore, from all methods 
presented here to characterize the morphology of the ionization 
fraction field, the power spectrum of the ionization fraction is the 
measure most closly related to observations. 

Shown in Figure [6] is the ionized fraction power spectrum, 
^xx{k), for the two simulations as a function of global ion- 



ization fraction on the abscissa. Instead of plotting A^(fc) = 
P{k) /{2iT^) against fc, we choose to plot it against 2.46/fc, the 
reason for this will become clear later. This colour contour plot 
thus shows the evolution of Axx{k) during reionization. For the 
53Mpc_g8.7_I30S model it can be seen that the power spectrum 
first peaks at scales of the order s = 2.46/fcniax ~ 0.8 Mpc. We 
expect this peak to be associated with the the size of ionized or 
neutral bubbles, for the following simple reason. On scales smaller 
than the bubbles, the correlation function £,xx{ri2) ~ {{x{ri) — 
Xv){x{r2) — Xy)) reduces to the constant value Xv{l — x^), while 
on scales much larger than the bubbles, the ionized fraction is 
uncorrelated, and the correlation function should approach zero 
jZaldarriaga et al..2004.) . This behavior for the correlation function 
implies that the power spectrum A^(fc) should approach zero at 
large and small scales, with a peak at the characteristic size of the 
bubbles. 

The first peak of Axx(fc) for a single spherical top-hat bubble 
of radius s would be located at fcmax ~ 2.46/s, which is why we 
use 2.46/fcniax to characterize the typical radius of regions. Indeed, 
comparison of the maxima of the spherical average model and the 
peaks of the power spectra shows that they are approximately re- 
lated by 2.5/fcmax ~ 4i?niax as can be seen by comparing Fig.|5] 
and Fig. [6] It should be noted that both, the spherical average and 
the power spectrum do not show a pronounced peak at all global 
ionization fractions, but are instead rather flat at higher global ion- 
ization fractions. This is in agreement with one of the trends that 
IZahn et al] ( 1201 Oh identify for all their tested reionization models. 
To better show this behaviour of the PS at higher global ionization 
fractions, we show the power spectrum of the fiducial simulation 
as a function of k for four different global ionization fractions as a 
line-plot in the left panel of Fig. |7] However, the other trend they 
identify, the shift towards smaller k with increasing lower ioniza- 
tion fraction is less apparent for our fiducial simulation. Instead, 
the form of the power spectrum of the fiducial simulation suggests 
that at global ionization fraction greater than (x) ~ 0.4 there are 
two main populations of bubbles: One at sizes where the PS ini- 
tially peaks, around 0.8 Mpc and the other with sizes around 6 Mpc, 
which is of the size of clusters of galaxies, i.e. the typical cluster- 
ing scale of the source halos at this epoch, not to be confused with 
the length scale which encompasses the mass of the higher-mass, 
virialized galaxy clusters familiar at lower redshift. The first scale 
results from suppression of sources after initial turn on of a source 
in a low-mass halo. Growth of the H II region is completely halted 
until a high-mass halo forms in that region. A less stringent sup- 
pression criterion would slow down the growth but probably not 
halt it entirely. The second scale results from merging of bubbles 
emerging from galaxies in the same cluster. These details in the 
size distribution are washed out in the spherical average distribu- 
tion. However, it should be noted that a scale of 6 Mpc is a sub- 
stantial fraction of a 53 Mpc box and therefore it is questionable if 
the sampling at this scale is high enough. At global ionization frac- 
tions larger than {x) ~ 0.4, it can be seen that there is considerable 
power on scales comparable to the box size. 

For the 53Mpc_gl.7_8.7S simulation it can be seen that there 
is more power on smaller scales but also that the power spectrum 
is flatter than the one of the fiducial model as there is no distinct 
peak below global ionization fractions of 20%. Further it can be 
seen that the slope of isochromatic lines is greater for this model 
than for the fiducial one. This means that the sizes of H II regions 
in the 53Mpc_gl.7_8.7S are growing faster with respect to global 
ionization fraction. This was also seen in the spherical average re- 
sults. The absence of the peak at smaller scales that is present in the 
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Figure 7. left panel: Power spectium of the fiducial simulation at several 
different global ionization fractions as indicated in the legend; right panel: 
Cross-correlation coefficient r^g as defined in Eq.|6]as a function of 2.46/fc 
at {x) = 0.05 and 0.5 (as indicated in the figure (by line thickness) for the 
fiducial simulation (gray lines) and 53Mpc_gl.7_8.7S (black lines). Note 
the lower correlation between ionized fraction and density for the fiducial 
simulation on all scales at low global ionization fractions. 

fiducial model, is due to the fact that already at a global ionization 
fraction of roughly 10%, the contribution from sources hosted by 
massive halos is about the same as the contribution from sources in 
low-mass halos, while this is true at about 70% global ionized frac- 
tion for our fiducial simulation. Therefore, the relative contribution 
from H II regions produced by sources in low mass halos in iso- 
lated cells is smaller. H II regions produced by sources in massive 
halos will grow continuously, explaining the flat distribution below 
several Mpc. 

The right panel of Fig. |7] shows the cross-correlation coeffi- 
cient of ionized fraction and density field, 

r^sik) = ^''[''^ ^.^ (6) 

at two different (x) ~ 0.05 and 0.5 for both simulations plotted 
against s = 2.46/fc. When r^s = (—1)1, the ionized fraction 
and density field are perfectly (anti-)correlated, while r^s — im- 
plies they are uncorrelated. As seen from the figure, the ionized 
fraction and density fields are nearly perfectly correlated on large 
scales, s > 8 Mpc. It can be seen that the scale s at which the 
correlation starts to decrease, is increasing with global ionization 
fraction. This is due to the fact that while the H II regions grow, 
they also start ionizing the voids. At very low global ionization 
fractions, represented here by (x) ~ 0.05, the correlation coeffi- 
cient for simulation 53Mpc_gl.7_8.7S is greater than the one for 
the fiducial simulation (especially at smaller s). This is expected 
because the ionizing radiation of the sources in the simulation with 
lower efficiencies can less easily "break out" of high density re- 
gions. Additionally, less efficient sources trace the high density re- 
gions better since clustered low mass sources are less suppressed: 
individually they form smaller H II regions and therefore do not 
suppress each other. At later stages of reionization this difference 
disappears: The simulation with low source efficiencies reaches the 
same global ionization fraction as the fiducial simulation at much 
later times when massive sources are more common. Those massive 
sources form bigger HII regions which also grow into the voids. 

3.3 Topology of Reionization 

Minkowski functionals have been used extensivel y in cosmology to 
characterize the to pology of large scale str ucture dGottetalJI 19861 : 
iMecke et al.lll994 ISchmalzing & Buchert"l997^ and also the non- 
Gaussianity of the cosmic microwave background (.Komatsu et al.l 
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I2OO9I) . Recent work has focused on using Minkowski functionals 
as a way to characterize the morphological structure of reionization 
dOleser et al. 2006; Lee et al. 2008). 

Both works focused on the topology of the H I density field. 
They showed the Euler Characteristic (or genus, respectively) as a 
function of neutral density for several different times (i.e. global 
ionization fractions) and concentrated on the increasing deviations 
from the typical curve of a gaussian random field. Here, we will 
take a complementary approach based upon the topology of the 
ionization fraction field, rather than the fluctuating neutral density 
field. Unlike the neutral density field, which takes values spread 
continuously over a very wide range, the ionized fraction field 
ranges only between and 1, and, ideally, there would essential 
be only two values possible to assign to any given point in space, 
either "neutral" values close to zero or "ionized" values close to 
unity. In that ideal case, the Euler Characteristic for the ionized 
fraction field would only be a function of time (or of the evolving 
globally-averaged ionized fraction) and be largely independent of 
the choice of ionization fraction threshold. 

W e follow the definition and notation of ISchmalzing et akl 
( Il996l) and Schinalzing & Buchert ( 1997). Consider a scalar func- 
tion /(x) defined at each point x G R^. The set Fth of all points 
X for which /(x) > /th defines bodies in three dimensional space. 
The zeroth Minkowski functional, Vq (/th), is simply the volume of 
those bodies: 

l^o(,fth) = / e[/th - ./■(x)]d^x, (7) 
Jv 

where O is the Heaviside step function. The next three Minkowski 
functionals are defined as surface integrals over the boundary of the 
bodies: 

V'i(/«0 = ^ / (x)d'A (8) 



V2{U) = 
^3(/th) = 




where Ri and R2 are the principal radii of curvature along the 
surface dFth- The first Minkowski functional is proportional to 
the integrated surface area. This and the zeroth Minkowski func- 
tional were used when calculating the size estimator 3V/A. The 
Minkowski functional V3, which is proportional to the integral of 
the Gaussian curvature over the surface, is also known as the Euler 
characteristic, and is equal to; 
Mparts — T^tunnels + ^cavities 

Applied to the ionization fraction field, two disconnected ionized 
cells would count as two parts, a ring-like ionized region constitutes 
a tunnel and one part and a neutral cell completely surround ed by 
ionized cells is a cavity. Table 1 in ISchmalzing et al.l l ll996l) gives 
an overview of different notations for the Minkowski functionals 
which only differ in constant factors. The better known quantity 
genus, g (nuinber of complete cuts one can make through the ob- 
ject without dividing it into disconnected parts), is related to the 
Euler characteristic by the simple relation g — 1~ VJj.lflFor exam- 
ple, a torus has V3 — 0, since it has zero total curvature, and has 

" This is true if one considers the Euler Characteristic of the volume de- 
fined by the set of points. Note that others consider the Euler characteristic 
of the surface of the set of points defining the volume. In this case, the 
Euler Characteristic X ^ factor of two greater, resulting in the relation 
to genus: x = 2(1 — g). This is consistent with the relation x{d^) = 
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Figure 8. Euler characteristic V3 as a function of global ionization frac- 
tion (x) and threshold value for the fiducial simulation (simulation 
53Mpc_gl.7_8.7) in the left (right) panel. The fast changes in the lower 
right comer are due to the implementation of photons that would travel dis- 
tances longer than the box-size, see explanation in the text. Note the high 
dependence on threshold value for simulation 38Mpc_gl.7_8.7. 

one part and one tunnel. A sphere, on the other hand, has V3 = 1, 
since it has one part and no tunnels, and positive total curvature. 

We oversample the ionization fraction fields before calculating 
the Euler characteristic. We do this to minimize critical connections 
of H I and H II regions. A critical connection is for example an ion- 
ized cell which is connected via an edge to another ionized cell in 
an otherwise neutral neighbourhood. Appendix|C]explains in more 
detail the problems that are involved. Important for the following 
analysis is to note that oversampling reduces ambiguities concern- 
ing the connectivity at a given threshold value but introduces higher 
dependences of V3 on the threshold value. 

In the left panel of Fig. [8] we show the evolution of the Euler 
characteristic V3 of the ionized fraction as a function of the thresh- 
old Xth for our (oversampled) fiducial simulation. We choose to 
normalize V3 by dividing by the box size, to have an easier com- 
parison when dealing with different box sizes. Between threshold 
values Xth = 0.2 and 0.6, the evolution of V3 is largely indepen- 
dent on the actual choice of a^th: V-g, rises to a maximum value at 
a mean ionization fraction of about 5%, after which V3 decreases 
and gets negative before 20% global ionization fraction is reached. 
It rises again after the ionization fraction passed 50% but never 
reaches positive values again. 

This behaviour can be qualitatively understood by consider- 
ing inside-out reionization in an approximately Gaussian density 
field. V3 for a Gaussian random field (or any monotone and steady 
function of it) as a f unction of threshold value, as plotted for ex- 
ample in figure 1 in ISchmalzing & BuchertI jl997h . shows in the 
second half a rise to positive values and decreases again. In Fig.[9l 
we show V3 of the density field from the 53 Mpc box simulation at 
redshift 2: ~ 26.1 as a function of density- threshold value Q The 

x(A)(l + (—1)''"^), where d is the dimension, A a d-dimensional body 
and d A its d — 1-dimensional surface, see equation (18) in Mecke et al. 
il994) . 

''^ V3 of the original density field shows an asymmetry between isolated re- 
gions and isolated cavities, for details see Appendix|c] A Gaussian smooth- 
ing with (T ~ 3 cells would be necessary to account for this. This removes 
all small scale structure, therefore V3 is substantially reduced. It also re- 
moves extreme over- and under- densities which is why the V3 -curve gets 
narrower Using sub-grid sampling enhances cells with intermediate densi- 
ties and therefore changes the distiibution away from gaussian distribution. 
Therefore we see deviations from the expected curve for a gaussian random 
field, which has | V3™ | / 1 ^3™="^ | = 1/(2 exp(-3/2)) ~ 2.2, as can b e 
calculated for example with equation (14) in lSchmalzing & Bucherll jl997l) . 




iog,„(p /<p>)| 



Figure 9. Euler Characteristic V3 of the density field as a function of density 
threshold value logjQ(p/ (p). Shown is V3 of the original density field at 
redshift 2 = 26.1 and V3 of oversampled fields as indicated in the figure. 
The inset shows V3 of the smoothed field with two different smoothing 
width, as indicated. 



fact that V3 of the ionization field does not show a rise to positive 
values followed by a decrease (see Fig.[8]or Fig.|C2), shows that if 
there exists a monotonic steady functional relation between density 
and time of ionization it does so only until a certain density: low- 
density areas (the voids in the density field) are ionized "before 
their time" and do not serve as positive contributions to Vz in the 
ionization field. This is another feature of inside-out reionization: 
the H II regions do eventually break out in the voids. 

It can be seen that V2. changes very fast at high global ion- 
ization fractions especially when choosing lower threshold values, 
see lower right comer in the left panel of Fig. [8] This is due to 
the way photons that would travel further than a box-distance are 
implemented in the simulations: The lost photons are collected and 
evenly distributed over all cells. Cells that only get ionized by those 
photons have small ionization fractions that only depend on their 
density. 

For very high values of Xth, the connectivity is reduced, there- 
fore V3 is greater due to the positive contribution from many more 
disconnected H II regions. This reduced connectivity is due to 
many partly ionized cells originating partly by the relatively small 
photon-output per halo per ionization time step compared to the av- 
erage number of atoms in a cell and partly by the oversampling that 
introduces additional partly ionized cells. 

It can be seen (Fig. [S] right panel) that the dependence of 
V3 on the threshold value is more pronounced for simulation 
53Mpc_gl.7_8.7S. From the EOF investigations (see Fig|4j, we 
know that this simulation has many cell-size H II regions. Also 
the simple estimates in Table 1 show that the smallest sources are 
not efficient enough to ionize their own cell, producing partly ion- 
ized cells. Therefore, a strong dependence on Sth is to be expected. 
This shows that the resolution for this choice of source efficiencies 
is not sufficient for doing topological investigations. Due to our 
choice of oversampling the data, choosing a higher threshold value 
(xth ~ 0.5) corresponds to underestimating the connection of H II 
regions. Choosing a lower threshold value might overestimate the 
connectivity. Comparing values of V3 at higher and lower threshold 
values, can be used as an indication of the sufficiency of resolution 
of the simulation. In the remainder of this work, we will mostly 
only show the evolution of V3 at threshold value xth = 0.5 of 
the oversampled data-fields and indicate cases where V3 is highly 
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dependent on the threshold value. Since the extrema of V3 of dif- 
ferent simulations can vary quite a bit, we choose to plot the func- 
tion /(V3) = RexA'T — Imv^Vs instead of just V3. In the inter- 
val [Vr^, V-rn ^ [f{Vm, f{Vmi / is bijective (i.e. the 
function has an inverse function), therefore we continue to refer to 
it as V3. 



4 BOX-SIZE AND RESOLUTION 

In this section, we investigate the effect of simulation volume size 
on the simulation and the effect of resolution on our analysis meth- 
ods. We compare our fiducial simulation to a simulation with the 
same source properties but in a bigger volume, 163Mpc_g8.7_130S. 
Before we do so, we analyze a smoothed version of the data of our 
fiducial simulation to test the effect of resolution on our analysis 
methods. We replace the ionization fraction data in each cell with 
the average over a three cell width volume centered on the cell in 
question, i.e. an average over 27 cells. We will refer to this as three- 
cell smoothing. This results in a resolution similar to the one in the 
163 Mpc simulation. It should be kept in mind that smoothing over 
three cells does not remove all structure smaller than three cells. 

In Sect. 3 we introduced three measures of size distribution 
and one estimate for the average bubble size. In Fig. [TO] we plot 
all four measures as a function of global ionization fraction for the 
smoothed version of 53Mpc_g8.7_130S. We also show curves of 
the peaks of the spherical average and the power spectrum. As a 
reference, these same curves (including the ?,V/A estimator) for the 
fiducial simulation are included as gray lines. We first concentrate 
on the iV/A estimate. As can be seen in the middle panel of FigfTOl 
(comparing the two solid lines), above a global ionization fraction 
(a;) ~ 0.1, the smoothed version yields larger values for ?,V/A. 
However, the difference is never greater than 20%. 

The effect on the spherical average distribution (same panel) is 
mainly a reduction of contribution from scales below 0.5 Mpc. Ad- 
ditionally, the peak of the spherical average (compare dot-dashed 
lines in the same panel) is slightly shifted towards larger scales for 
the smoothed data. 

A somewhat contrary effect can be seen in the FOP size dis- 
tribution: The contribution from scales below s ~ 0.3 Mpc is en- 
hanced in the smoothed data. This can be explained as follows. 
If a larger ionized structure is elongated (i.e. no structure in 2 di- 
mensions) and very inhomogeneous in its ionization fraction, then 
3-dimensional smoothing would break up the structure in smaller 
parts. 

In the right panel of Fig.[TO]we show the power spectrum of 
the three-cell smoothed data, its maximum curve and the maxima 
from the fiducial simulation without smoothing. It can be seen, sim- 
ilar to the spherical average distribution, that power on scales be- 
low a ~ 0.5 is removed. Also, it can be seen, that at higher global 
ionization fraction, the distinct peak at scales around 0.8 Mpc di- 
minishes while the peak at s ~ 6 Mpc is as pronounced as in the 
unsmoothed data. 

Since smoothing reduces the small scales and therefore re- 
duces the critical (vertex/edge) H II/H I region-connections, over- 
sampling the smoothed data does not change V3 more than 10%. 
It can be seen that the form of the evolution of V3 for the fiducial 
simulation stays roughly the same even for an eleven-cell smooth- 
ing of he data. However, since the smaller scales are smoothed out, 
the total amplitude of V3 is reduced. 

Equipped with an idea about which effects can be due to 
the changed resolution, we now turn to the larger volume sim- 
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Figure 11. Evolution of V3 for smoothed versions (degree of smoothing 
as indicated in the figure by line-thickness) of the fiducial simulation at 
threshold value I'th = 0.5. Note the similarity of the behavior of the curves 
for different smoothing lengths. 



ulation. In Fig. [T2] we plot all 4 size measures for simulation 
I63Mpc_g8.7_130S. We concentrate first on the ?)V/A size esti- 
mate, see the middle panel of that figure (black solid line) and com- 
pare it to the iV/A estimate of the fiducial simulation (gray solid 
line): except for very low global ionization fractions i{x) < 0.08) 
where the estimate for the fiducial simulation is of the order of the 
cell size of the 163 Mpc simulation volume, the two curves almost 
coincide up to a global ionization fraction (x) ~ 0.3, after which 
the scale in the 163 Mpc simulation grows faster. 

The SPA distribution shows a similar behaviour as 3V/A: be- 
low ionization fractions (x) ~ 0.08, the SPA peaks of the fiducial 
simulation are at scales comparable to the cell size of the 163 Mpc 
simulation while the peaks of 163 Mpc simulation are slightly 
larger (compare the dot-dashed lines in the same panel). The fidu- 
cial simulation shows also a wider distribution with contributions 
from smaller as well as from larger scales, see left panel in Fig|5] 
Between (a;) ~ 0.08 — 0.3 the evolution of SPA distribution of 
the two simulation is very similar. At larger global ionization frac- 
tions, the scale of the SPA peak is growing faster in the 163 Mpc 
simulation. 

The power spectrum also shows that, below global ionization 
fractions (z) ~ 0.3, despite the different resolution, the simula- 
tions with different simulation volumes agree well. At global ion- 
ization fractions (x) > 0.4 there is considerable power on scales 
that are not captured in the 53 Mpc simulation which might explain 
the shift found by the other size distribution estimates. The power 
spectrum suggests that the 163 Mpc simulation captures the most 
relevant scales involved in reionization (that is H II region sizes up 
to roughly (x) ~ 0.5, H I region sizes above this) since there is 
little power on scales above s ~ 30 Mpc. A noticeable difference 
between the simulations is the lack of the peak at scales around 
0.8 Mpc in the large simulation volume. As found earlier by the 
smoothing test, this may be a resolution effect. The peak at scales 
s ~ 6 Mpc is not as clear in the 163 Mpc volume as in the 53 Mpc 
volume. After (x) ~ 0.5 it shifts to larger scales. 

The FOF size distribution (left panel of Fig.[T2]and left panel 
of Fig|4l, of the two simulations look on first sight very different. 
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Figure 10. Color plots of evolution of size distributions for the three-cell smoothed data of simulation 53Mpc_g8.7_130S. The left panel shows the FOF 
distribution (color coded as in Fig.|4). The middle panel shows the SPA (color coded as in Fig.|5) and its peaks (black dot-dashed line) together with the ZV/A 
estimate (black solid line). The right panel shows the PS (color coded as in Fig.|6) and its peaks (black dot-dashed line). For comparison, the gray lines show 
the corresponding measures of the fiducial simulation. Indicated on the ordinates are the cell and simulation volume size for the 53 Mpc simulation (>) and 
the cell size of the 163 Mpc simulation volume (-). 
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Figure 12. Color plots of evolution of size distributions for simulation 163Mpc_g8.7_130S, from left to right: FOF, SPA (including ZV/A, solid black line) 
and PS. Color coding and line-styles have the same meaning as in Fig. 1101 



It can be easily understood why: the smallest scale of H II regions 
in the small simulation volume is smaller than that in the large sim- 
ulation volume. Looking at the cell- volume limits which are indi- 
cated on the abscissa, it can be seen that the additional population 
of small scale H II regions present in the 53 Mpc simulation is be- 
low the cell-size of the 163 Mpc simulation and therefore below 
its resolution limit. All those H II regions are partly ionized in the 
large simulation volume. Therefore, some of them which are more 
ionized than a;th = 0.5 appear as additional population in the FOF 
distribution of the large simulation volume at scales of its cell size. 
At 10% global ionization fraction, it can be seen that the 53 Mpc 
simulation has slightly larger large-scales than the 163 Mpc sim- 
ulation, s ~ 4Mpc and s ~ 5Mpc, respectively. Partly ionized 
cells that are ionized below the threshold value in the larger simu- 
lation volume and therefore do not count as belonging to the H II 
region cannot account completely for this difference in size. This 
slight mismatch in size between large and small simulation volume 
might be an artifact from our implementation of suppression: since 
we suppress all sources inside a cell that has a mass averaged ion- 
ization fraction Xm larger than 10%, the volume in which sources 
are suppressed can be overestimated in simulations with larger cell 
sizes. 



While the size distribution in the small simulation volume 
shows a gap after (a::) ~ 0.1, the gap emerges first at {x) ~ 0.25 
in the large simulation volume. This is because a single region in 
any of the bins that are empty in the small simulation volume but 
populated in the large one, would already exceed the contribution it 
makes in the small volume. This sampling effect was already men- 
tioned in the previous section. 

At 30 % global ionization fraction, the ionized volume of the 
single large connected region in the 163 Mpc simulation is approx- 
imately 10% of the total simulation volume which is larger than the 
size of the 53 Mpc volume. Similarly, the volume of the largest con- 
nected ionized region in the 53 Mpc simulation is also about 10% 
of the total volume. The fact that the largest connected region is a 
constant fraction of the simulation volume already at 30% ionized 
fraction suggests that this region pervades the whole simulation 
volume. This statement is strengthened by the Euler Characteristic 
of these simulations: V3 is already highly negative at (x) ~ 0.3 (for 
the 163 Mpc simulation this is only true for lower threshold values), 
see left panel of Fig.[T3] It should be noted that Vi of the 163 Mpc 
simulation shows very similar evolution to V3 of the 53 Mpc sim- 
ulation (for a low threshold value for the 163 Mpc volume). The 
fact that V3 in the large simulation volume is highly dependent on 
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Figure 13. Left panel: Evolution of Euler characteristic V3 at tliresliold 
value Xth = 0.5 and 0.3 as indicated in the figure (by line style) for both 
g8.7_130S simulations, gray lines indicating the 53 Mpc simulation and 
black lines the 163 Mpc simulation. 

Right panel: Cross congelation of the density field with ionization fraction 
of both g8.7_130S simulations (with the same color coding as in the left 
panel) at two different global ionization fractions as indicated in the figure 
(by line style). 

threshold value, shows that the resolution is not sufficient to use V3 
as a reliable analysis tool: the ambiguity as to whether regions are 
connected or not is too high, as can be seen comparing the curves 
for different threshold values in the left panel of Fig[T3] However, 
it can be seen that the effect of including "lost" photons, is much 
smaller in the large simulation volume: relatively fewer photons 
travel distances longer than a box distance in the large simulation 
volume. 

The cross correlation between ionized fraction and density at 
low global ionization fractions is almost identical for both simula- 
tion volumes as the thin lines in the right panel of Fig. [T3] show. 
It differs more at higher global ionization fractions where it shifts 
towards larger scales for the 163 Mpc simulation, probably indi- 
cating that the scales dominating the ionization field at that global 
ionization fraction exceed the size of the 53 Mpc simulation. 



5 PHYSICAL PARAMETERS 

Our efficiency parameter is a product of the efficiency of star for- 
mation, production of ionizing photons per stellar atom (related to 
the initial stellar mass function) and the escape fraction of the pho- 
tons from the galactic halo into the intergalactic medium. All these 
quantities are not very well constrained at present. Also the effi- 
ciency of suppression due to Jeans-mass filtering can be different 
from the simple on-of f function as implemented in our simulations 
with suppression (c.f [ McOuinn et alj |2007| : iMesinger & Diikstral 
I2OO8I : lOkamoto et alj|2008f) . To studv the effect of our simplified 
suppression model, we consider in this section some extreme sce- 
narios and use the methods described above to investigate the effect 
on the scales and the topology of the emerging H II regions. 

5.1 Minimum mass of halos hosting sources witli escaping 
ionizing radiation 

In this subsection we investigate how a change in the source pop- 
ulation affects the simulation. In simulations 53Mpc_uvS_le9 and 
53Mpc_gl0.4_0 only halos more massive than IQ^ Mq host sources 
that emit ionizing radiation into the IGM. The former simulation is 
constructed such that the number of released ionizing photons in 
every time step (after the formation of the first massive halos) is 
the same as for our fiducial simulation. The sum of all photons that 



were released in the fiducial simulation by low mass halos up to 
the time at which the first massive halo is forming, is emitted ad- 
ditionally in the first time step after which the first massive source 
has formed. Since the number of forming halos increases exponen- 
tially, the fraction of additionally released photons in this first time 
step is only about half of the total released photons at that time step. 
As pointed out in section |2l the resulting source efficiency is vari- 
able with time, shown in Fig[T] This also means that the minimum 
number of photons released by one source during one time step is 
decreasing with increasing (x). Therefore, also the minimum size 
for H II regions decreases to {x) ~ 0.25. This can be seen most 
clearly in the FOF size distribution (Fig. [141 IsfO and in the power 
spectrum (Fig. [14] right). To avoid this effect we performed simu- 
lation 53Mpc_gl0.4_0 which has a different ionization history from 
our fiducial simulation, but the source efficiency of the high mass 
sources is chosen such that overlap occurs at roughly the same time, 
as can be seen in table 1. 

The FOF size distribution shows that individual H II regions 
grow larger before merging with the largest H II region in both, the 
53Mpc_uvS_le9 and 53Mpc_fl0.4_0 simulations than in the fidu- 
cial one; the gap in the distribution is smaller This is due to the 
greater average distance between high mass sources. The space 
in between the large H II regions is neutral, without any ionized 
spots. Therefore, each individual H II region can grow bigger be- 
fore merging. 

Below global ionization fractions (x) ~ 0.2, the evolution of 
sizes in the 53Mpc_uvS_le9 simulation is dominated by the first 
H II regions emerging around the highly efficient first sources. At 
higher global ionization fractions the size-evolution is very simi- 
lar to simulation 53Mpc_gl0.4_0. Therefore, we concentrate in the 
following on the 53Mpc_gl0.4_0 simulation. 

Compared to our fiducial simulation, the 3V/A estimate of 
simulation 53Mpc_gl0.4_0 (middle panel of Fig[T5j suggests an 
average bubble scale about a factor 3 greater at all global ionization 
fractions we consider here. Also the spherical average distribution 
clearly shows this shift to larger scales. This is best visible when 
comparing the peak-scales (dot-dashed curves in the same panel). 
The power spectrum (Fig. [15] right panel), reveals that it is only 
a shift to bigger scales below global ionization fractions of about 
20%. The H II regions that form first seem to be larger than the 
ones in the fiducial simulation. Later, it is rather a lack of small 
scales; notably the peak at scales s ~ 0.8 Mpc is absent. This is 
not surprising as we identified the peak to be due to the suppres- 
sion of sources in low mass halos. The suppresion is responsible 
for halting the growth of the H II regions formed by these sources 
completely. The peak at scales s ~ 6 Mpc is still there. However, 
there is more power on scales that are not captured by the 53 Mpc 
simulation volume. 

The Euler Characteristic for simulation 53Mpc_gl0.4_0, see 
light gray lines in the left panel of Fig. [TS] is in total much flatter 
than the Euler Characteristic for the fiducial simulation. Since we 
saw in the FOF distribution and the power spectrum that there is 
not much contribution from scales that are hardly resolved, this is 
very unlikely to be due to unresolved scales. We therefore conclude 
that the bigger H II regions around rare sources result in a less 
complex topology with fewer tunnels and cavities. We saw above 
that feeding back diffuse photons into the volume, affects V3 at 
high global ionization fractions. The fraction of those photons for 
simulation 53Mpc_gl0.4_0 is already 0.2 at 80% global ionization 
fraction. Therefore, the evolution of Vs beyond (x) ~ 0.8 might be 
dominated by the effect of these photons, as described in section[3] 
It can be seen that V3 for the low threshold value is at some points 
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Figure 14. Color plots of evolution of size distributions for simulation 53Mpc_uvS_le9, from left to right: FOF, SPA (including 3V/A, solid black line) and 
PS. Color coding and line-styles have the same meaning as in Fig.EO] 




Figure 15. Color plots of evolution of size distributions for simulation 53Mpc_gl0.4_0, from left to right: FOF, SPA (including 3V/A, solid black line) and 
PS. Color coding and line-styles have the same meaning as in Fig. 1101 



considerably different from the value at xth = 0.5. However, since 
we plot the square root of V3, the differences at values close to 
are amplified. The actual difference is never greater than 20% of 
the maximum value. 



5.2 Source suppression vs. low efficiency 

To test the effect of source suppression in regions where the IGM 
is ionized, we compare our fiducial simulation (with instantaneous 
complete suppression of sources in low mass halos in ionized re- 
gions) to two simulations without suppression: 53Mpc_g8.7_130 
which has the same source efficiencies but which ends much ear- 
lier due to the many more released photons and 53Mpc_g0.4_5.3 
which has substantially lower source efficiencies to end at roughly 
the same time as the fiducial simulation. 

The comparison of the FOF size distributions between model 
53Mpc.g8.7_130 and 53Mpc.g0.4_5.3 (see left panels in Fig. [16] 
and Fig. I17t shows that the model with the lower source efficien- 
cies shows more very small H II regions than the fiducial simu- 
lation, similar to 53Mpc_gl.7_8.7S. The simulation with the same 
source efficiencies as the fiducial simulation but without suppres- 
sion shows less small H II regions than the fiducial simulation, be- 
cause each individual source forming is active longer and so con- 
tinuously grows its H II region. Also, clustered sources in the same 
or neighboring cells support the growth of their joined H II region. 



The spherical average distribution, see middle panels of Fig. 
[76] and Fig. [17] for the higher and lower efficiency simulations 
without suppression, respectively, show a very similar evolution. 
The rate at which the average size grows seems to be the same, 
but the scale is shifted towards larger scales for the simulation 
53Mpc_g8.7_I30 by about a factor 1.5. Also, the 3V/A size esti- 
mates suggests an almost constant shift to larger scales. 

The power spectra, see right panels of the same figures show 
more power on small scales (below s ~ 0.5) and less power on 
large scales (above s ~ 6) for 53Mpc_g0.4_5.3 than for simula- 
tion 53Mpc_g8.7_130. However, up to global ionization fraction 
(x) ~ 0.5, the peak at scales around s ~ 6 Mpc is present in 
both simulations. The peak at scales s ~ 0.8 Mpc is absent in both 
simulations. 

It should be noted that the size distributions found by all three 
methods as well as the 3V/A size estimator of the 53Mpc_g8.7_I30 
simulation are very similar to the ones from 53Mpc_gI0.4_0. How- 
ever, there is a very significant shift in time between these two sim- 
ulations. The reason for their similarity if compared at equal (x) 
might be that most halos which first reach masses above 10** Mq, 
reach masses above 10^ Mq at accordingly lower redshifts. There- 
fore, the same halos hosting sources which turn on as sources in low 
mass halos in the former simulation, turn on as sources in more 
massive halos at a later time. This means that if suppression of 
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Figure 16. Color plots of evolution of size distributions for simulation 53Mpc_g8.7_130, from left to right: FOF, SPA (including ZV/A, solid black line) and 
PS. Color coding and line-styles have the same meaning as in Fig. 1 101 




Figure 17. Color plots of evolution of size distributions for simulation 53Mpc_g0.4_5.3, from left to right; FOF, SPA (including 3V/A, sohd black line) and 
PS. Color coding and line-styles have the same meaning as in Fig. 1 101 



sources is not important, the sources "shaping" reionization are the 
ones in halos with the lowest mass that can form luminous sources. 

The Euler Characteristic (left panel in Fig. I18t shows that 
Va never becomes negative for both, simulation 53Mpc_g0.4_5.3 
(dark gray lines) and simulation 53Mpc_g8.7_130 (black lines): 
There are many more disconnected H II regions than neutral tun- 
nels through the ionized regions. For a lower threshold value, 
xth = 0.3, Vs goes mildly negative. The high dependence of sim- 
ulation 53Mpc_g8.7_I30 on threshold value is somewhat surpris- 
ing since the smallest H II regions are larger than for the fiducial 
simulation. However, the total number of H II and H I regions is 
small (only a few hundred compared to a few thousand for the fidu- 
cial simulation), therefore a few critical connections are enough to 
introduce a strong dependence of V3 on threshold value. For the 
53Mpc_g0.4_5.3 simulation it can further be seen that the maxi- 
mum V3 is about a factor of four larger (at low threshold value 
which means most probably overestimating the connections of re- 
gions) than for the fiducial simulation. There are more discon- 
nected H II regions due to their smaller minimum sizes. Simula- 
tion 53Mpc_g8.7_130 reaches the same global ionization fraction 
as simulation 53Mpc_g0.4_5.3 at a much higher redshift, at a time 
where fewer halos have formed. The maximum V3 is smaller for 
the 53Mpc_g8.7_130 simulation because the individual H II regions 
grow bigger and merge earlier in terms of global ionization fraction. 
However, the big differences between V3 at different threshold val- 
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Figure 18. Left panel: Va of simulation 53Mpc.gl0.4.0, 53Mpc.g8.7.130 
and 53Mpc_g0.4_5.3, as indicated in the figure (by line color) at threshold 
values a;tij = 0.5 and 0.3, as indicated in the figure (by line style). 
Right panel: Cross-correlation coefficient for the same simulations at dif- 
ferent (z), as indicated in the figure (by line style). 



ues show that these simulations do not have sufficient resolution, 
therefore the ambiguity of connections is too high. 

Fig. [TS] (right panel) shows that the correlation between ion- 
ized fraction and overdensity at small global ionization fractions, 
here represented by (a;) ~ 0.05, is flatter for simulations without 
suppression, or simulations without any suppressible sources, like 
simulation 53Mpc_gl0.4_0: on large scales, the correlation is lower 
than for simulations with suppression like our fiducial simulation; 



© 2010 RAS, MNRAS 001.[TH20l 



16 Friedrich et al. 



on small scales, relatively larger. This behaviour can be interpreted 
together with the results from the Euler Characteristic in the follow- 
ing way: The lack of suppression leads to earlier break outs of the 
ionizing radiation into the low density regions between low mass 
sources which reduces the number of neutral tunnels and lowers the 
cross-correlation on larger scales. At the same time, not suppress- 
ing sources in partly ionized regions leads to complete ionization of 
that region which increases the cross-correlation on smaller scales. 

At higher global ionization fractions, here represented by 
{x) ~ 0.5, these differences disappear because almost all sources 
in low mass halos are suppressed in the simulations with suppres- 
sion and the reionization process is dominated by the sources in 
high mass halos. 

iMcOuinn et alj jioOTh tested the effect of suppression for the 
case of equally efficient high and low- mass sources. However, their 
simulations do not resolve halos below l(fyiQ. They use an ana- 
lytic prescription to include unresolved halos above the H I atomic 
cooling mass. Among other things they found that even their most 
drastic (instantaneous) suppression model (complete suppression 
of sources in low mass halos) yields an ionization field with sim- 
ilar morphology to a simulation without suppression (their simu- 
lations F3 and SI, respectively). Additionally they tested a sim- 
ulation with higher efficiencies for low mass sources but without 
suppression (their model S2). To realize this they chose a mass de- 
pendent source efficiency factor. Very roughly, this translates into 
(lOVlO^)"^/^ ^ 5 times more efficient low (10^ > AZ/Mq > 
10*) than high (M/M© > 10^) mass sources in terms of our step- 
function efficiency assignment method (our fiducial model has 15 
times more efficient sources in low mass halos). They found that the 
effect of boosting the efficiency of low mass sources on ionization 
field morphology is rather small. From these two tests (F-set vs. Si 
and Si vs. S2) they concluded that suppression, even in the case of 
high efficient low mass sources, does not affect the morphology of 
the ionization field. We can confirm that the ionization field mor- 
phology does not change much if the efficiency of low mass sources 
is boosted, as long as the total photon output of the least luminous 
sources that form the smallest HII regions is held roughly con- 
stant (compare for example 53Mpc_gl0.4_0 and 53Mpc_g8.7_130, 
Fig.[T5]and Fig.[T6] respecively). However, comparing our fiducial 
simulation (with suppression) to simulations 53Mpc_g8.7_130 and 
53Mpc_g0.4_5.3 (both without suppression), we find large differ- 
ences in the topology and in the average size distributions of ion- 
ized regions, as discussed above. This is most probably due to the 
early break out of H II regions formed around sources in low-mass 
halos in the case of no suppression and high low-mass halo source 
efficiency. 



6 CONCLUSIONS 

We have explored different methods for characterizing the scales 
and topology of complex ionization fraction fields produced by 
simulations of cosmic reionization. For characterizing the length 
scales or sizes of HII regions, we used three methods that give a dis- 
tribution of scales; the FOF method, the spherical average method 
and the power spectrum of the ionization fraction field. In addition 
we proposed a single valued measure of average size of HII region, 
given by the ratio of the volume to the surface of all regions. For 
characterizing the topology we employed the Euler Characteristic 
or third Minkowski functional, Vs of the ionization fraction field. 

The nature of the size distribution of H II regions can be 
viewed to be a matter of definition. Applying a literal definition 



leads to the FOF approach, in which the HII regions are considered 
to be connected regions of space. Because the topology of reion- 
ization can be quite complex, as seen from the Euler characteristic, 
this definition does not lend itself easily to analytical modeling, and 
connecting the FOF size distribution to other scale estimators, such 
as the power spectrum, is by no means trivial. For the FOF method, 
what is lost in its complexity is gained in the detailed description 
of the reionization process that it provides. Although most of the 
volume is already at quite low global ionization fractions contained 
in one large connected region, there is a wealth of information con- 
tained in the number and sizes of the smaller bubbles, which only 
occupy a small fraction of the volume. As we found in section]?] the 
FOF size distribution is affected by resolution in terms of cell size 
and one has to take this into account when interpreting the results 
and comparing simulations with different resolution. In principle, 
the FOF method yields the maximum size for isolated ionized re- 
gions, before they start to merge into the largest one. However, in 
practice, this scale is dependent on the sampling of the distribution 
function and therefore on the size of the simulation volume. 

The spherical average method, on the other hand, gives dis- 
tributions which are much smoother and can more easily be con- 
nected to analytical models jZahn et al .12001) . The results are more 
sensitive to the large scale H II regions which constitute the main 
contribution to the global ionization fraction. Due to its averaging 
nature, it washes out the details. Consider, for example, a collection 
of two types of spheres with scales si = a and S2 = & > a; the 
spherical average method will only reveal the two distinct scales if 

< 1/3. Also, the spherical average alone tends to substantially 
underestimate the sizes of H II regions, as shown in Fig. lAII bv our 
toy model. For using it as a size estimator on par with other meth- 
ods, the spherical average scales should be multiplied by a factor 
4. 

When the universe is mostly neutral, the peak of the ionization 
power spectrum is related to the size distribution of ionized regions. 
For a single top-hat sphere, the first peak in the power spectrum 
is related to its radius by fcmax = 2.46/r and we choose to use 
2.46/fc as a size estimator when comparing to the results of other 
methods. The power spectrum produces size distributions roughly 
comparable to those from the spherical average method. However, 
the advantage of the power spectrum over the spherical average 
method is that it does not wash out the details of the distribution. As 
mentioned earlier, of all size estimators discussed in this paper, the 
power spectrum is the one most related to upcoming observations 
as o ne component in the e xpansion of the 21 cm power spectrum 
(e.g. lFurlanetto et al.ll2006l) . 

The ratio of the zeroth and first Minkowski functionals can be 
used to define a mean radius of HII regions, the iV/A estimate. 
This estimate is a surface weighted average. It generally gives re- 
sults consistent with the maximum of the spherical average, and 
when it falls below that there is a large fraction of small bubbles 
(which dominate the surface). 

The Euler characteristic V-j, of the ionization fraction field of- 
fers a rich description of the evolution of the topology of reion- 
ization. Taking our fiducial simulation as an example, in the early 
stages of reionization the value of V3 is positive as the topology 
is dominated by a large number of isolated H II regions. However, 
already beyond global ionization fractions of roughly 20%, V3 be- 
comes highly negative, indicating a complex topology of connected 
HII regions with tunnels. This is consistent with the distribution 
from the FOF method, which shows that already at {x) ~ 0.3 the 
main contribution to ionized fraction comes from only one large 
connected region which pervades most of the simulation volume. 
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As the ionization fronts around stellar sources are quite thin, 
ideally V3 should not depend much on the chosen threshold value, 
and the most interesting aspect is the change of 1/3 with time. In 
this sense topology studies of HII regions differ from those of den- 
sity fields, where the variation of Vi with threshold value is used 
to characterize the field at a given time. However, in practice the 
resolution of the ionization field may not be sufficient to achieve 
sharp fronts, and partially ionized cells will occur. The result is dif- 
ferent values of V3 for different threshold values and even differ- 
ent evolutions of V3 at different threshold values. More seriously 
is the interaction of this effect with the definition of connectiv- 
ity/adjacency used when calculating the Euler Characteristic. We 
proposed a new test for establishing how robust the derived values 
of V3 are to a change of adjacency. In this one compares the an- 
swer with the value obtained for a field and threshold value which 
have their sign inverted. Using this we showed that sub-sampling 
or smoothing is generally required to obtain consistent results, and 
that for fields with a large fraction of partially ionized cells, it can 
be difficult to get consistent results. 

We subsequently applied the size estimators and Euler Char- 
acteristic to study differences and similarities between different 
reionization simulations. Comparing identical simulations in two 
different volume sizes, 163 Mpc and 53 Mpc, shows that below 
global ionization fractions of 30% the average scales of H II re- 
gions are roughly the same for both simulations. Beyond that the 
size distributions in the larger volume start to contain scales be- 
yond those available in the smaller one. Another manifastation of 
this is that the largest connected region found by the FOF-method 
for both simulation volume sizes is 10% of the each of the simula- 
tion volumes already at 30% global ionization fraction. Therefore, 
this largest H II region is about a factor 30 larger in the large sim- 
ulation volume: A region of this size does not fit into the 53 Mpc 
simulation volume. 

Even earlier there are differences between the two EOF dis- 
tributions showing an absence of H II regions with volumes of a 
few hundred Mpc"^ in the small box. This is because relatively iso- 
lated overdensity regions (surrounded by larger voids) are missing 
in the small box due to numerical variance. Since the contribution 
from those scales is very small, the effect on the average scales is 
negligible. In general, simple size estimates may be biased due to 
missing scales (e.g. intermediate scales missing in the small box 
due to the smaller sampling volume, small scales missing in big 
box due to resolution). If one is not interested in very small scales, 
the lower resolution in the big box appears not to be a problem. The 
fact that the position of the peak of the power-spectrum at all times 
(global ionization fractions) is well below the box-size scale for the 
163 Mpc box, indicates that no larger simulation boxes are needed 
to follow the evolution of the peak position. 

Comparing simulations with and without suppression (either 
with only high mass sources, or with sources in low mass halos not 
suppressed in ionized regions) shows that the ones without suppres- 
sion typically have a much less complex ionization fraction field 
topology and a more steady growth of average H II region sizes. We 
thus find significant differences between the cases with and with- 
out suppression. Eor the simulations with sources only in high mass 
halos, this can be explained as follows: the individual H II regions 
can grow much bigger before merging, due to the larger average 
distances between high-mass sources. Eor the simulation without 
suppression and high efficiency sources a similar argument holds, 
but the whole process takes place at much earlier redshifts. Eor the 
simulation without suppression and with low efficiency sources the 
explanation is again similar, but now applied to smaller scales in- 



stead of at earlier times. This similarity between the cases without 
suppression is due to the statistical nature of the density field. If 
suppression is not important, the sources in halos with the lowest 
mass still capable of forming sources are the ones shaping reion- 
ization. 

By imposing an external photon-budget on an independently 
evolved source population, one can in principle separate the effect 
of the source population from that of the reionization history. How- 
ever, by necessity this implies an evolution of the source efficien- 
cies, which in the case we studied seriously affected the evolution 
sizes of HII regions. Up to a global ionization fraction (x) ~ 0.2 
the size distribution is dominated by the size of H II regions gener- 
ated by the first generation of sources. In the later stages of reion- 
ization, the morphology of the ionization field is very similar to a 
simulation with a fixed source efficiency for high mass sources. 

As outlined in the introduction, we concentrated on the anal- 
ysis of the ionization fraction fields and its evolution. Applying 
the various analysis methods to the future observations of the red- 
shifted 2 1 cm signal is in principle possible, but requires sufficient 
sensitivity to image the signal at different frequencies. The first 
generation of telescopes is not expected to be able to do this, but 
the planned Square Kilometer Array (SKA)EI should be. Compared 
to the simulation results, the observations will have the additional 
complications of noise and limited spatial resolution (below even 
that of our 163 Mpc simulation). As we have shown, resolution 
effects should be treated with care, especially for the Euler Charac- 
teristic, and noise peaks can obviously also bias the topology deter- 
mination. Still, characterizing the morphology of HII regions in the 
data will be important as they trace the mass and thus the emerging 
cosmic web. We leave the application of the various size/scale and 
topology estimates to mock observational data to a future paper. 
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APPENDIX A: SPHERICAL AVERAGE SIZE 
DISTRIBUTION OF A LOG NORMAL DISTRIBUTION 

The spherical average method was described bv lZahn et al.l ( l2007h . 
In this technique, each cell in the computational volume is con- 
sidered to be in an ionized region if a sphere centered on that cell 
has a mean ionized fraction greater than a given threshold, usually 
xth = 0.9. The size of the H II region to which it belongs is taken 
to be the largest such sphere for which the condition is met. 

In the following, we assume that gas is either fully ionized 
or fully neutral, and that all ionized bubbles are non-overlapping 
spheres with a volume-weighted distribution dP/dR, so that 
P{R + dR) — P{R) is the fraction of the ionized volume that 
lies within bubbles with radii between R and R + dR. What bubble 
distribution, dPsm/dR, would be obtained by using the spherical 
average method? To simplify further, we will take the threshold for 
the spherical average, a;th, to be arbitrarily close to unity, so that a 
point is considered to be within an ionized sphere of a given radius 
only if all the matter in that sphere is ionized. For a single ionized 
sphere of radius r, dPsm{R)/dR is the surface of a sphere with 
radius r ~ R, normalised by the integral of the surface of spheres 



with radii from to r: dPsin{R) /dR = 3(r - 
this as W{r, R), then 



dPs, 



dR 



drW{r, R) 



dP 

dr ' 



'7?)7r'^.Ifwe define 



(Al) 



is the bubble size distribution obtained by the spherical average 
method for the real distribution dP/dR. The lower limit of the in- 
tegral is R because only spheres which are larger than R can con- 
tribute to the spherical average bubble distribution at R (because 
xth — > 1); the largest ionized sphere that can be drawn around any 
given point is always smaller than or equal in radius to the actual 
ionized sphere in which it lies. 

Shown in Fig. lAll are rdP/dr and the corresponding 
RdPsm/dR for three log-normal distributions of bubble sizes 

^2- 



dP 
din r 



1 



: exp - 



(ln(r)-ln((r))- 



2a2 



(A2) 



with different a. As can be seen from the figure, the spherical av- 
erage tends to change the true bubble distribution in two ways. 
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Figure Al. Spherical average method applied to three log-normal distri- 
bution of non-overlapping spherical H II regions. The right (soUd) curves 
are the actual log-normal distributions given by equation jA2t for different 
values of cr, while the left (dashed) curves show the result which would be 
obtained on the same distribution using the spherical average method. As 
CT — > and dP/dR approaches a delta function, dPsm/dR approaches 
the kernel function W{R, {R)) 

First, it smooths the actual bubble distribution with the kernel 
W{r,R). Second, it lowers the value of the mean bubble radius, 
Rnv = J RdP/dR. In our simple toy model, the mean bubble size 
obtained by the spherical average method is always 1/4 of the ac- 
tual one. Our toy model is admittedly crude, most notably in the 
assumption of a threshold Sth = 1 and spherical H II regions. A 
lower value of xth would allow small pockets of neutral gas to be 
attributed to large ionized regions. In fact, for the case where x = 1 
and 2: = in ionized and neutral regions, respectively, lower val- 
ues of a;th lead to an overestimate of the volume which is ionized, 
leading to a violation of the normalization condition, 

/ diR-TB = (A3) 
Jo dR 

In most cases this overestimate is not very large. A lower value of 
Xth would also yield larger H II regions, but this effect has been 
shown to be rather modest jZahn et al]|2007h . The assumption of 
spherical symmetry is a conservative one, however, in the sense that 
it provides a lower limit to how much the spherical average method 
underestimates the "true" H II regions sizes. This is because the 
spherical average method is sensitive to the smallest dimension of 
the region in which it lies: the radius is larger than the smallest 
dimension, the part of the sphere lying in that direction would lie 
outside the region, and the average ionized fraction would no longer 
be above the threshold for the region to be considered ionized. 



APPENDIX B: ADDITIONAL SIZE MEASURE 

iMesinger & Furlanettd j2007h used another technique to charac- 
terize sizes of ionized (and neutral) regions in their binary ioniza- 
tion fields: They chose a large number of random points; for each 
point they check if it is ionized (neutral); from each ionized (neu- 
tral) point they chose a random direction and measure the distance 



from this point to the nearest ionized-neural (neutral-ionized) tran- 
sition boundary along that line of sight. In the following, we use a 
method similar to the one in Mesinger & Furlanetto (2007). How- 
ever, instead of choosing a random direction, we, for simplicity, 
only choose between the 3 principle axis in one of the two possi- 
ble orientations. The differences in size distribution obtained with 
this method only matters for very large and rare H II regions. Small 
regions should be abundant enough that the orientation should not 
play a dominant role. Since we have to deal with continuos ion- 
ization fields, we have to introduce two parameters: Above which 
ionization fraction is a random point ionized? What is the limit for 
the transition boundary along the line of sight? In the following, 
we consider a point as ionized if its ionization fraction is above 
Kth; we count each point along the line of sight ionized as long 
as its ionization fraction is above a^um. In Fig lBll we plot the size 
distribution curves obtained with this method for the two simula- 
tions we use in Section [3] 53Mpc_(78.7_130S' (our fiducial sim- 
ulation, solid lines) and SSMpc.gl.T.S.TS (dashed lines) at two 
different global ionization fractions, {x) ~ 0.1 (thin lines) and 
(a;) ~ 0.4 (thick lines). Qualitatively, we find the same result as 
with the other size measures: Initially the size distribution peaks 
earlier for 53A/pc_(;1.7_8.75; at higher (x), the typical size seems 
to be larger in this simulation than for the fiducial simulation. 

With respect to the peak position of the (uncorrected) SPA, 
especially at low global ionization fraction, the peak of this size 
distribution is shifted towards larger scales. This is consistent with 
the findings in appendix A and the fact that the size distribution of 
a single sphere, found with this method, would peak at the radius of 
the sphere; therefore, this method yields in theory better results than 
the SPA. However, the disadvantage is the dependence on two pa- 
rameters if the field of investigation is not a binary field. We tested 
the effect of these parameters (compare different colours in Fig lBlb 
and find large variations for different parameter combinations. 



APPENDIX C: ON THE CALCULATION OF THE EULER 
CHARACTERISTIC 

The estimation of V-j, of a three dimensional field sampled on a fi- 
nite set of grid points is dependent on the chosen adjacency-pair, if 
the structure is of the same order as the cell sizes of the grid. An ad- 
jacency pair is for example (26,6) which means that cells above the 
threshold (foreground cells) have 26 neighbours and cells below the 
threshold (background cells) have 6 neighbours. Because of this de- 
pendence on the choice of adjacency pair, we choose to oversample 
the data to avoid differences in the treatment of isolated/connected 
H II and isolated/connected H I regions. This introduces a stronger 
dependency on the choice of the threshold value. 

For calculating the Euler Characteristic V-j, we use part 
of a program developed b y T. Buchert and J. Schmalzing 
dSchmalzing & Bucherl|[T997h . The algorithm we use counts the 
vertices {V), edges {E), faces {F) and lattice cells (C) of the 
foreground cells and calculat es V3 according to equation (10) in 
ISchmalzing & Bucherll l\wk : 

Vi = V - E + F -C (CI) 

In terms of adjacencies, this is equivalent to assigning 26 neigh- 
bor cells to any foreground cell and 6 neighbor cells to any back- 
ground cell. This adjacency pair ensures the 3 dimensional Jordan 
curve theorem (which basically means that an edge-connection of 
foreground cells cannot at the same time be a connection for back- 
ground cells). .Ohser et aL t200Z) showed that this adjacency pair 
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Figure Bl. Size distribution curves for simulations g8.7_130S' (solid lines) 
an d (71.7-8. 75 (da shed lines) using a method similar to the one described 
in iMesinger & Furlanettol ^20071) with three different combination of pa- 
rameters Xtij and xiiui (line colour) at two different global ionization frac- 
tions (hne thickness). Also indicated are the maxima of the curves and the 
maxima of the SPA for the same simulations at the same global ionization 
fractions. Note the large effect of the choice of parameters x^^ and xw^^ 
on the position of the curve maximum. Note also the shift of the maxima 
of the SPA towards smaller scales. The maximum of SPA for simulation 
gl.7_8.7S' at global ionization fraction {x) ~ 0.1 is at 0.1 Mpc and there- 
fore not visible in this plot. 

is complementary which means that V3 of the background is the 
same as that of the foreground. If the structure in the data-cube has 
contributions smaller or of the same size than is sampled by the 
grid-cells, structure of lower dimensions than 3 can arise. This can 
cause inconsistencies in the approximation of V3 of the set sampled 
at the grid-points resulting in a violation of the complementarity. 
This is for example visible in the {V-i, 5th) plot in Fig.|9] The first 
peak is mainly due to disconnected under-dense regions which are 
below the density threshold value (5th. If two of those regions are 
connected via an edge or a vertex, they count as two disconnected 
"cavities" since those cells are background cells and have only 6 
neighbours. The second peak is somewhat smaller than the first 
one (although theoretically the peaks should be equal). This peak 
is mostly due to disconnected over-dense regions that are above the 
threshold value (5th and therefore have 26 neighbours each. Over- 
densities that are connected via an edge or a vertex count as one 
connect ed region , there fore their contribution to V3 is smaller. 

Ohs eretai] ( I2OO2I) showed that the bias of the approximation 
depends on the choice of adjacencies. Inverting the data cube of 
ionization fraction, so that at every grid point, x' = —x, com- 
puting Vi at x[i^ = -xth and comparing Vi{x[y,) to V3(xth), is 
equivalent to changing the adjacencies from (26/6) to (6/26). 

For example, consider a 3'' cube CI, with Cl(l,l,2) = 
Cl(2, 2, 2) — 1 and everywhere else and a cube C2 = —CI. 
Then, y3(Cl,xth) / 14(C2,-a;th): V3(Cl,a:th ^ I) ^ U - 
23 + 12 — 2 = 1, where I G (0, 1) (round brackets denoting open 
intervals) while V3(C2,a;th = -I) = (i9)-(C) + (C)-(v-2) = 2, 
where ■£), ^, ( and are the numbers of vertices, edges, faces and 
lattice cells of the cube C(l : 3, 1 : 3, 1 : 3) = 0, the sum of which 
is since periodicity is assumed. 

Therefore, we use V3{—x, — a;th) as a check for how good 
the approximation of V3 is. To minimize the asymmetry due to the 
chosen adjacency, we oversample the data. For the example of the 
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Figure CI. Two-dimensional example to demonstrate the effect of over- 
sampling the data: original structure (left panel) compared to the result from 
oversampling (right panel). While the original structure counts as one con- 
nected region independent of threshold value (z^h G (0, 1)), the structure 
in the right panel counts as two disconnected regions if Xth > 0.5, and 
as one connected region otherwise. If the entries would be inversed (multi- 
plied by -1), the structure in the left panel would count as two disconnected 
cavities while the structure in the right panel would have the same depen- 
dence on the threshold value as before: the cavities count as disconnected if 
l^thl > 0.5, and as one connected cavity otherwise. 

cubes CI and C2 from above, this means that V3(C1, aJth ~ I) ~ 
V3{C2,xn, = -0 = 1 for ' G (0,0.5) and V3(Cl,Xth = I) = 
V3{C2,Xth = -0 = 2 for / G [0.5,1). Obviously this results 
in a higher dependency on the chosen threshold value, as this sim- 
ple example demonstrates, see also Fig. lCll for a two-dimensional 
example. 

As a second example, we present here the Euler Characteris- 
tic of the fiducial simulation and simulation 53Mpc_gl.7_8.7S to 
demonstrate the effect of oversampling and the choice of threshold 
value. As explained above, calculating V3(— a;th) of the inverted 
field (the data field multiplied by —1, in the following indicated by 
a "-"), is equivalent to changing foreground to background cells and 
vice versa. This is the same as changing the adjacencies. Therefore, 
comparing V3 (+, a;th) to V3 (— , — xth) is a test for the dependence 
on adjacencies. 

Fig |C2| demonstrates the dependence on the chosen adjacency 
pair when analyzing the data-fields without smoothing or oversam- 
pling (compare V3(+,a;th = 0.5) to V3(— ,a;th = —0.5) for 
both simulations) and the lower dependence on the chosen adja- 
cencies after oversampling the data (compare V3"^(+,a;th = 0.5) 
to V3°^(— ,2:th = —0.5)). It also demonstrates the relatively low 
dependence on threshold value a::th for the fiducial simulation 
and the high dependence on Xth for 53Mpc_gl.7_8.7S (compare 
Vr{+, xtb = 0.5) t oT/3°''(+.Xth = 0.3)). 

Others (e.g. Gleser et al I2OO6I) chose to smooth the data with 
a gaussian kernel to remove lower-dimensional parts. This has the 
disadvantage that some of the small scale structures are suppressed. 
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Figure C2. Evolution of Euler Characteristic V3 for the original field of the 
fiducial simulation (lower gray thick line) and simulation 53Mpc_gl.7_8.7S 
(lower black thick line) at threshold value Xth = 0.5. The upper thick 
lines are the corresponding Euler Characteristic of the inversed fields at 
threshold value Xth = —0.5. The thin lines correspond to the oversam- 
pled fields of the fiducial simulation (lower gray thin line) and simulation 
53Mpc_gl.7_8.7S (lower black thin line) at threshold value x^^^ = 0.5. The 
upper thin lines are the corresonding Euler Characteristic of the inversed 
oversampled fields at threshold value x^^^ = —0.5. The gray (fiducial sim- 
ulation) and black (53Mpc_gl.7_8.7S) lines with crosses indicate the Euler 
Characteristic of the oversampled fields at threshold value x^^ = 0.3. Note 
the good agreement for V^^{xt\i = 0.5) and V^^{xx^ = —0.5) and the 
low dependence on threshold value for the fiducial simulation. 
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